home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d26 / matlab.arc / PDE < prev    next >
Encoding:
Text File  |  1990-09-23  |  1.2 KB  |  41 lines
  1. //Conductivity example.
  2. //Parameters
  3.   rho   //radius of cylindrical inclusion
  4.   n     //number of terms in solution
  5.   m     //number of boundary points
  6. //initialize operation counts
  7.   flop = <0,0>;
  8. //initialize variables
  9.   m1 = round(m/3);
  10.   m2 = m - m1;
  11.   pi = 4.0*atan(1.0);
  12. //generate points in Cartesian coordinates
  13.   //right hand edge
  14.   for i = 1:m1, x(i) = 1; y(i) = (1-rho)*(i-1)/(m1-1);
  15.   //top edge
  16.   for i = m2+1:m, x(i) = (1-rho)*(m-i)/(m-m2-1); y(i) = 1;
  17.   //circular edge
  18.   for i = m1+1:m2, t = (pi/2)*(i-m1)/(m2-m1+1); ...
  19.     x(i) = 1-rho*sin(t); y(i) = 1-rho*cos(t);
  20. //convert to polar coordinates
  21.   for i = 1:m-1, th(i) = atan(y(i)/x(i)); ...
  22.     r(i) = sqrt(x(i)**2+y(i)**2);
  23.   th(m) = pi/2; r(m) = 1;
  24. //generate matrix
  25.   //Dirichlet conditions
  26.     for i = 1:m2, for j = 1:n, k = 2*j-1; ...
  27.       a(i,j) = r(i)**k*cos(k*th(i));
  28.   //Neumann conditions
  29.     for i = m2+1:m, for j = 1:n, k = 2*j-1; ...
  30.       a(i,j) = k*r(i)**(k-1)*sin((k-1)*th(i));
  31. //generate right hand side
  32.   for i = 1:m2, b(i) = 1;
  33.   for i = m2+1:m, b(i) = 0;
  34. //solve for coefficients
  35.   c = A\b
  36. //compute effective conductivity
  37.   c(2:2:n) = -c(2:2:n);
  38.   sigma = sum(c)
  39. //output total operation count
  40.   ops = flop(2)
  41.